load motor.mat;
XtestNorm = motorNormalize(Xtest);
XtrainNorm = motorNormalize(Xtrain);

meanPredictionsPoly = zeros(100,1);
alpha = logspace(-8,0,100);
alpha = alpha';
lambda = zeros(size(alpha), 1);
for a = 1:size(lambda)
   lambda(a, 1) = 400*alpha(a, 1);
end
RMSEpolynomialFor100 = zeros(size(alpha), 3);
RMSEradialFor100 = zeros(size(alpha), 3);
meanPredictionPolynomial = zeros(size(alpha), 3);
meanPredictionRadial = zeros(size(alpha), 3);

for a = 1:size(alpha)
    RMSEpolynomialFor100(a, 1) = alpha(a, 1);
    RMSEradialFor100(a, 1) = alpha(a, 1);
    meanPredictionPolynomial(a, 1) = alpha(a, 1);
    meanPredictionRadial(a, 1) = alpha(a, 1);
end



%sigma^2/lambda = 1/alpha
%sigma^2 * alpha = lambda
%400*alpha = lambda

numDataPoints = size(XtrainNorm,1);
for lambdaIndex=1:size(alpha)
    Mindex = 100;
    powerIndex = Mindex;
    
    designMatrix = ones(numDataPoints,powerIndex+1);
    for rowDesign=1:numDataPoints
        for colDesign=2:powerIndex+1
            designMatrix(rowDesign,colDesign)=XtrainNorm(rowDesign,1).^(colDesign-1);
        end
    end
    w = (designMatrix' * designMatrix* (lambda(lambdaIndex,1)*eye(Mindex+1)))\designMatrix' * Ytrain;
    
    YPredict = zeros(size(XtestNorm,1),1);
    for xtestIndex=1:size(XtestNorm,1)
        YPredict(xtestIndex,1) =  w(1,1);
        for xPower=2:powerIndex+1
            YPredict(xtestIndex,1) = YPredict(xtestIndex,1) + XtestNorm(xtestIndex,1).^(xPower-1)*w(xPower,1);
        end
        RMSEpolynomialFor100(lambdaIndex,2) = RMSEpolynomialFor100(lambdaIndex,2) + (Ytest(xtestIndex) - YPredict(xtestIndex,1)).^2;
    end
    RMSEpolynomialFor100(lambdaIndex,2) = sqrt(RMSEpolynomialFor100(lambdaIndex,2) / size(XtestNorm, 1));
    meanPredictionPolynomial(lambdaIndex,2) = sum(YPredict)/size(XtestNorm,1);
    
    YPredict = zeros(size(XtrainNorm,1),1);
    for xtrainIndex=1:size(XtrainNorm,1)
        YPredict(xtrainIndex,1) =  w(1,1);
        for xPower=2:powerIndex+1
            YPredict(xtrainIndex,1) = YPredict(xtrainIndex,1) + XtrainNorm(xtrainIndex,1).^(xPower-1)*w(xPower,1);
        end
        RMSEpolynomialFor100(lambdaIndex,3) = RMSEpolynomialFor100(lambdaIndex,3) + (Ytrain(xtrainIndex) - YPredict(xtrainIndex,1)).^2;
    end
    RMSEpolynomialFor100(lambdaIndex,3) = sqrt(RMSEpolynomialFor100(lambdaIndex,3) / size(XtrainNorm, 1));
    meanPredictionPolynomial(lambdaIndex,3) = sum(YPredict)/size(XtrainNorm,1);

end

meanPredictionsRadial = zeros(100,1);

numDataPoints = size(XtrainNorm,1);
for lambdaIndex=1:size(alpha)
    Mindex = 100;
    mu = linspace(-1,1,Mindex);
    sigma = 3*(mu(2)-mu(1));
    
    designMatrix = ones(numDataPoints,powerIndex+1);
    for rowDesign=1:numDataPoints
        for colDesign=2:powerIndex+1
            designMatrix(rowDesign,colDesign)=exp(-(XtrainNorm(rowDesign,1)-mu(colDesign-1))^2 / 2*sigma^2);
        end
    end
    w = (designMatrix' * designMatrix* (lambda(lambdaIndex,1)*eye(Mindex+1)))\designMatrix' * Ytrain;

    testPhi = ones(size(XtestNorm,1),Mindex+1);
    for x=1:size(XtestNorm,1)
        for j = 2:Mindex+1
            testPhi(x,j) = exp(-(XtestNorm(x)-mu(j-1))^2 / 2*sigma^2);
        end
    end
    
    YPredict = zeros(size(XtestNorm,1),1);
    for xtestIndex=1:size(XtestNorm,1)
        YPredict(xtestIndex,1) = testPhi(xtestIndex,:)*w;
        RMSEradialFor100(lambdaIndex,2) = RMSEradialFor100(lambdaIndex,2) + (Ytest(xtestIndex) - YPredict(xtestIndex,1)).^2;
    end
    RMSEradialFor100(lambdaIndex,2) = sqrt(RMSEradialFor100(lambdaIndex,2) / size(XtestNorm, 1));
    meanPredictionRadial(lambdaIndex,2) = sum(YPredict)/size(XtestNorm,1);
    
    trainPhi = ones(size(XtrainNorm,1),Mindex+1);
    for x=1:size(XtrainNorm,1)
        for j = 2:Mindex+1
            trainPhi(x,j) = exp(-(XtrainNorm(x)-mu(j-1))^2 / 2*sigma^2);
        end
    end
    
    YPredict = zeros(size(XtrainNorm,1),1);
    for xtrainIndex=1:size(XtrainNorm,1)
        YPredict(xtrainIndex,1) = trainPhi(xtrainIndex,:)*w;
        RMSEradialFor100(lambdaIndex,3) = RMSEradialFor100(lambdaIndex,3) + (Ytrain(xtrainIndex) - YPredict(xtrainIndex,1)).^2;
    end
    RMSEradialFor100(lambdaIndex,3) = sqrt(RMSEradialFor100(lambdaIndex,3) / size(XtrainNorm, 1));
    meanPredictionRadial(lambdaIndex,3) = sum(YPredict)/size(XtrainNorm,1);

end

figure(271)
semilogx(RMSEpolynomialFor100(:,1),RMSEpolynomialFor100(:,3))
figure(272)
semilogx(RMSEpolynomialFor100(:,1),RMSEpolynomialFor100(:,2))
figure(273)
semilogx(RMSEradialFor100(:,1),RMSEradialFor100(:,3))
figure(274)
semilogx(RMSEradialFor100(:,1),RMSEradialFor100(:,2))


